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ABSTRACT 

In 2008 the blazar Markarian 421 entered a very active phase and was one of the brightest sources in the sky 
at TeV energies, showing frequent flaring episodes. Using the data of ARGO-YBJ, a full coverage air shower 
detector located at Yangbajing (4300 m a.s.l., Tibet, China), we monitored the source at gamma ray energies E 
> 0.3 TeV during the whole year. The observed flux was variable, with the strongest flares in March and June, 
in correlation with X-ray enhanced activity. While during specific episodes the TeV flux could be several times 
larger than the Crab Nebula one, the average emission from day 41 to 180 was almost twice the Crab level, 
with an integral flux of (3.6±0.6) x 10~ u photons cm -2 s -1 for energies E > 1 TeV, and decreased afterwards. 

This paper concentrates on the flares occurred in the first half of June. This period has been deeply studied 
from optical to 100 MeV gamma rays, and partially up to TeV energies, since the moonlight hampered the 
Cherenkov telescope observations during the most intense part of the emission. Our data complete these obser- 
vations, with the detection of a signal with a statistical significance of 3.8 standard deviations on June 11-13, 
corresponding to a gamma ray flux about 6 times larger than the Crab one above 1 TeV. 

The reconstructed differential spectrum, corrected for the intergalactic absorption, can be represented by 
a power law with an index a = — 2.1^oj extending up to several TeV. The spectrum slope is fully consistent 
with previous observations reporting a correlation between the flux and the spectral index, suggesting that this 
property is maintained in different epochs and characterizes the source emission processes. 
Subject headings: BL Lacertae objects: individual (Markarian 421) - gamma rays: observations 
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1. INTRODUCTION 

The blazar Markarian 421 has been the first extragalac- 
tic source observ ed at gamma ray energy E>500 GeV 
dPunchetal. 19921) . It belongs to the radio-loud Active 
Galactic Nuclei (AGNs) sub class of BL Lacertae objects, 
characterized by a non thermal spectrum extending up to high 
energies and by a rapid flux variability at nearly all wave- 
lenghts. 

To date about 30 BL Lacs have been detected at very high 
energies (VHE, E>100 GeV) and Mrk421 is the closest one 
(z = 0.031). Its relatively small distance makes it one of the 
best studied TeV gamma ray sources. Since its discovery this 
object played a significant role in the discussion concerning 
both the emission processes in AGNs and the attenuation of 
TeV gamma rays in the extragalactic space. 

It is now widely recognised that the BL Lac radiation origi- 
nates in a relativistic jet pointing at a small angle to the line of 
sight and that it is amplified by relativistic effects, explaining 
both the strong high energy emission and its rapid variability. 

Usually the BL Lacs energy density spectra have two 
broad band components, the first one peaking in the in- 
frared to X-ray region, the second one in th e Me V- TeV 
range dSambruna et al. 19961: iFossati et al. 19981) . Mrk421 is 
classified as a High-energy peaked BL Lac (HBL), show- 
ing the peaks in the X-ray and VHE regions, respectively 
( Pado vani & Giommi 19951) . 

The low energy component is commonly believed to orig- 
inate as synchrotron emission from relativistic electrons gy- 
rating in the magnetic field of the jet plasma, while the ori- 
gin of the second one is still unclear. Many models pro- 
pose that gamma rays are produced in Inverse Compton 
scattering of synchrotron (Synchrotron Self-Compton, SSC) 
or ambient photons (external Compton, EC) off the same 
electron population that produces the syn chrotron radiation 
dGhisellini et al. 19981: iDermeretal. 19921) . Alternatively, in 
the "hadronic" models, gamma rays are emitted as sychrotron 
radiation of extremely energetic protons, or by secondary par- 
ticle s produced by proto ns interacting with some target mate- 
rial dMiicke et al. 20031) . 

Flaring activity of Mrk421 at VHE energies has been ob- 
served with variability time scales ranging from minutes 
to months, and many multiwavelength campaigns have re- 
vealed a strong correlation of gamma rays with X-rays, 
that can be easily i nterpreted in terms of the SSC model 
dFossati et al. 20081: I Wagner 2008T) . In addition some data 
have shown significant variations of the TeV spectrum slope 
during different activity phases, and an evident correla- 
tion between the spe ctral hardness and the flux intensity 
dKrennrich et al. 20021) . 

The simultaneous observation at different wavelengths is of 
great importance since may provide unique information about 
the source properties and the radiation processes. 

A set of measurements dDonn arumma et al. 2009) covering 
12 decades of energy, from optical to TeV gamma rays, was 
performed during the strong flaring activity in the first half of 
June 2008 by different detectors: WEBT (optical R-band), 
UVOT (UV), RXTE/ASM (soft X-rays), SWIFT (soft and 
hard X-rays), AGILE (hard X-rays and gamma rays) and the 
Cherenkov telescopes VERITAS and MAGIC (VHE gamma 
rays). These data allowed for a deep analysis of the broad 
band energy spectrum as well as for the study of time correla- 
tions among the fluxes in different energy ranges. 

In this period two flaring episodes were reported, the first 



one on June 4-8, observed from optical to TeV gamma rays, 
the second one, larger and harder, on June 10-14, observed 
from optical to 100 MeV gamma rays. Using the multifre- 
quency data, Donnarumma et al. (2009) derived the Spectral 
Energy Distribution (SED) for June 6, that shows the typical 
two humps shape. In the framework of the SSC model, ac- 
cording to the authors, the second hump intensity, that reached 
a flux of about 3 x 10 _1 1 photons cm" 2 s" 1 at energies E >400 
GeV (i.e. about 3.5 times the Crab Nebula emission in the 
same energy range) seems to indicate that the variability is due 
to the hardening/softening of the electron spectrum, and not to 
the increase/decrease of the electron density. Their model pre- 
dicts for the second flare the Inverse Compton hump slightly 
shifted towards higher energies and a VHE flux a factor >2 
larger with respect to the first one. 

Unfortunately there were no VHE data included in their 
multiwavelength analysis after June 8 because the moonlight 
hampered the Cherenkov telescopes measurements. The VHE 
observation was actually made for such a very important flar- 
ing episode by the ARGO-YBJ experiment. 

The ARGO-YBJ experiment, located at the Yangbajing 
Cosmic Ray Laboratory (Tibet, P.R. China, 4300 m a.s.L, 30° 
06' 38" N, 90° 31' 50" E), since December 2007 is perform- 
ing a continuous monitoring of the sky in the declination band 
from -10° to 70°. 

In this paper we present our observation of Mrk421 in flar- 
ing state during 2008. After a summary of the data collected 
during the most active phase (February-June), we focus our 
discussion on the results obtained for the June flares in the 
framework of the Donnarumma et al. (2009) findings. 

2. THE ARGO-YBJ EXPERIMENT 

The ARGO-YBJ detector is constituted by a central car- 
pet ~74x 78 m 2 , made of a single layer of Resistive Plate 
Chambers (RPCs) with ^92% of active area, sorrounded by 
a partially instrumented (^20%) area up to ~100xll0 m 2 . 
The apparatus has a modular structure, the basic data ac- 
quisition element being a cluster (5.7x7.6 m 2 ), made of 12 
RPCs (2.8x1.25 m 2 ). Each chamber is read by 80 strips of 
6.75x61.8 cm 2 (the spatial pixels), logically organized in 10 
independent pads of 55.6x61.8 cm 2 which are individually 
acquired and represent the time pixels of the detector. The 
full detector i s made of 153 clu sters for a total active surface 
of ~6600 m 2 dAielli et al. 20061) . 

ARGO-YBJ operates in two independent acquisi- 
tion modes: th e shower mode and the scaler mode 
dAielli et al. 20081) . In the following we refer to the 
data recorded in shower mode. A simple, yet powerful, 
electron ic logic has been i mplemented to build an inclusive 
trigger dAloisio et al. 20041) . This logic is based on a time 
correlation between the pad signals depending on their 
relative distance. In this way, all the shower events giving 
a number of fired pads N pal j > N„, g in the central carpet in 
a time window of 420 ns generate the trigger. This trigger 
can work with high efficiency down to N,„ s =20, keeping 
negligible the rate of random coincidences. 

The time of each fired pad in a window of 2 /isec around the 
trigger time and its location are recorded and used to recon- 
struct the position of the shower core and the arr ival direction 
of the primary particle dDi Sciascio et al. 20071) . In order to 
perform the time calibratio n of the 18360 pad s, a software 
method has been developed dAielli et al. 2009aT) . 

The detector is in stable data taking with the trigger con- 
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dition N,„g=20 and a duty cycle > 85%. The trigger rate is 
^3.6 kHz with a dead time of 4%. 

3. DETECTOR PERFORMANCE 

The angular resolution and the pointing accuracy of the de- 
tector have been evaluated by using the Moon shadow, i.e. 
the deficit of cosmic rays in the Moon direction. The shape 
of the shadow provides a measurement of the detector Point 
Spread Function (PSF), and its position allows the individu- 
ation of possible pointing biases. ARGO-YBJ observes the 
Moon shadow with a sensitivity of about 10 standard devia- 
tions per month for events with a multiplicity N pa d >40 and 
zenith angle 9 <50°, corresponding to a proton median energy 
E p ~1 .8 TeV dDi Sciascio et al. 20081) . 

According to the Moon shadow data, the PSF of the de- 
tector is Gaussian for N pa d >100, while for lower multiplic- 
ities it can be described with an additional Gaussian, which 
contributes for about 20%. When the PSF is a Gaussian with 
r.m.s. a, the opening angle ip containing ~71.5% of the events 
maximizes the signal to background ratio for a point source 
with a uniform background, and is equal to 1.58 a. 

The semi-aperture ip is found to be 2.59°± 0.16°, 1.30°± 
0.14° and 1.04°± 0.12° for N pad >40, 100 and 300, respec- 
tively, in agreement with expectations from Monte Carlo sim- 
ulations. 

This measured angular resolution refers to cosmic ray- 
induced air showers. The angular resolution for 7-induced 
events has been evaluated by simulations and results smaller 
by ^30-40%, depending on N pa d, due to the better defined 
time profile of the showers. 

The relation between the observed pad/strip multiplicity 
spectrum and the primary energy spectrum has been studied 
with cosmic ray showers, by means of a f ull Monte Carlo s im- 
ulation, including the CORSIKA code dHecketal. 19981) to 
describe the shower development i n the atmosphe re, and a 
code based on the GEANT package (1GEANT 19931) to simu- 
late the detector response. Primary particles have been sam- 
pled from the energy spectra of Horandel (2003). The mea- 
sured strip multiplicity spectrum is in good agreemen t with 
the one predicted by the simulation (lAielli et al. 2009bT) . 

4. DATA ANALYSIS 

The dataset for the analysis of Mrk42 1 presented in this pa- 
per contains all showers with N pai j >40 and zenith angle less 
than 40°. No event selection based on the shower core posi- 
tion and no gamma-hadron discrimination have been applied 
in this work. 

A sky map in celestial coordinates (right ascension and dec- 
lination) with 0.1°x0.1° bin size, centered on the source lo- 
cation, is filled with the detected events. In order to extract 
the excess of 7 rays coming from the source, the cosmic ray 
background must be carefully estimated and subtracted from 
the event map. 

The b ackground is evaluated with the time swapping 
method dAlexandreas et al. 1993b . For each detected event, 
N "fake" events are generated by replacing the original ar- 
rival time with new ones, randomly selected from a buffer that 
spans a time T of data taking. We chose T ~ 3 hours to min- 
imize the systematic effects due to the environmental param- 
eters variations. Changing the time, the fake events maintain 
the same declination of the original event, but have a differ- 
ent right ascension. With these events a new sky map (back- 
ground map) is built. The number of fake events generated for 
each event is N = T(hr)x 15 x cos(<5)/2V>, where ip is the radius 



of the observational window in degrees (see below) and S is 
the declination of the source. In this way the average number 
of fake events falling in the observational window is ~1. 

The two maps are then "integrated" over a circular area of 
radius ip, i.e. every bin is filled with the content of all bins 
whose center has an angular distance less than ip from its cen- 
ter, with ip = 1.7°, 0.9° and 0.6° for N pad >40, 100 and 300, 
respectively. 

Finally the integrated background map is subtracted to the 
corresponding integrated event map, obtaining the "source 
map", where for every bin the statistical significance of the 
excess is calculated. 

With this procedure, however, since also the source events 
are used in the time swapping procedure, the obtained back- 
ground at the source position is slightly overestimated, and the 
signal underestimated. This underestimation increases with 
the observational window size, ranging from ~4 to 10% of the 
signal, depending on the N pai j interval. The observed event 
rate is then corrected using the appropriate factor. 

The whole procedure has been tested with the Crab Nebula, 
the standard candle for VHE astronomy. At the Yangbajing 
latitude the Crab culminates at a zenith angle 9 a ,i m = 8.1° and 
it is observable every day for 5.8 hours with a zenith angle 
(9<40°. 

The Crab Nebula was observed from 2007 December 13 
to 2009 August 8, for a total of 3150 on-source hours, ob- 
taining a signal with a statistical significance of 7.6 standard 
deviations for N pa d >40. The average number of gamma rays 
detected per day is 156.6±20.6 for N pad >40. 

To evaluate the energy spectrum, we simulate a source in 
the sky following the daily path of the Crab Nebula, and es- 
timate the number of events expected in different N p „d inter- 
vals, as a function of the spectrum parameters. 

Assuming a power law spectrum in the 0.1-80 TeV en- 
ergy range, the best fit to the data is dN/dE = (4.1 ±0.6) 
x 10" 11 (E/l TeV)- 2 ' 7±0 ' 2 photons cnr V 1 TeV" 1 , in agree- 
ment with our previous measurement dVernetto et al. 2009) 
and with observations by other detectors, such as H.E.S.S. 
dAharonian et al. 20061). MAGIC ([Albert et al. 20081) and Ti- 
bet AS-^ dAmenomori et al. 20091) . This result confirms the 
reliability of the simulation procedure and of the energy cali- 
bration of the detector. 

Concerning the energy range sampled in the Crab Nebula 
measurement, about 84% of the detected events comes from 
primary photons of energies greater than 300 GeV, while only 
8% comes from primaries above 10 TeV. 

The same analysis was performed for Mrk421. This source 
culminates at the ARGO-YBJ location at a zenith angle 
Ocuim = 8.1°, and it is observable every day for 6.4 hours with 
a zenith angle 9 < 40°. 

We evaluate the Mrk421 spectrum from day 41 to 180 of 
2008, when the X-ray flux showed the most intense flares. In 
this period (754 observation hours) the signal has a statistical 
significance of 5.8 standard deviations. 

We assume a power law spectrum multiplied by an expo- 
nential factor e~ r<£) to take into account the absorption of 
gamma rays in the Extragalactic Background Light (EBL), 
with the values of the optical depth t(E) given by Raue & 
Mazin(2008). 

The best fit spectrum obtained is: dN/dE = (3.0±0.5) x 
10- n (E/1.5 TeV)- 2 ' 4±0 - 3 e" r(£) photons cm" 2 s" 1 TeV" 1 . The 
integral flux above 1 TeV is (3.6±0.6) x 10~ n photons cm" 2 
s , almost twice the Crab Nebula one, i.e. 2.1 x 10~ n pho- 
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tons cm" 2 s" 1 , according to Aharonian et al. (2006). 

The values of the spectral index and of the gamma ray flux 
averaged over this 140 days period, support the correlation be- 
tween spectral hardness and flux intensity reported by Kren- 
nrich et al. (2002) and Albert et al. (2007), based on observa- 
tions of Mrk421 in different activity states. 

A complete account of the observations on Mrk421 will be 
reported in a dedicated paper. 

5. THE JUNE 2008 FLARES 

As mentioned in the Introduction, two different flares have 
been observed from Mrk421 in June 2008, the first one peak- 
ing in X-rays on June 4-6 and the second one on June 11-13. 
Concerning VHE gamma rays, Cherenkov telescopes data are 
available only for the first flare. An energy spectrum for E 
>400 GeV has been provided by VERITAS for June 6. 

Since the ARGO-YBJ sensitivity does not allow the obser- 
vation of a flux a few times larger than the Crab one in only 
one day (i.e. during one transit of the source in the detector 
field of view), we integrated the measurement over 3 days. 

Fig. 1 shows the rate of events with N pat j > 100 observed by 
ARGO-YBJ from June 3 to June 15, averaged over 3 days, 
compared with the X-ray flux measured by RXTE/ASJVQ in 
the 2-10 keV energy range. A correlation between the gamma 
ray and X-ray light curves is clearly visible. During the days 
June 11-13, when the maximum of the second flare occurred, 
the excess of events from Mrk421 reached a statistical signif- 
icance of 3.8 standard deviations. 

Beside the statistical error, this measurement could be af- 
fected by a systematic uncertainty due to the background eval- 
uation, that is the most delicate step of the analysis. In order 
to estimate this effect, a completely different procedure for 
the background calculation has bee n implemented, using th e 
so-called equi-zenith angle method ( Amenomori et al. 2005). 
In this method the events collected in 10 off-source windows 
of the same size of the on-source window, and aligned on 
both sides of the same zenith angle belt, are used to obtain the 
background. A detailed study of the two methods in the same 
sky region has shown that on average they give significances 
of the excesses consistent within 0.7 standard deviations, cor- 
responding to about 20% uncertainty on the flux estimate of 
the observed signal. 

The event rate as a function of the minimum pad multi- 
plicity, obtained integrating the data during June 4-6 for the 
first flare (17.9 hours) and during June 11-13 for the second 
one (18.2 hours), is shown in Fig. 2. On the same plot, the 
two solid lines represent the expected rates according to the 
Donnarumma et al. (2009) model, obtanined by a simulation 
procedure. The SED proposed by this model has been cor- 
rected for the EBL absorption using the parameters given by 
Raue & Mazin (2008) in order to have the flux at Earth. Then, 
using the absorbed spectrum, we simulated a source moving 
along the Mrk421 path on the sky, and evaluated the number 
of events expected in the detector, for different N pai i thresh- 
olds. The complete simulation procedure (which includes the 
gamma ray showers propagation in the atmosphere, and the 
detector response) has been tested evaluating the Crab Nebula 
flux, as shown in the previous section. 

In the limit of the statistical accuracy of this measurement, 
our data suggest for both flares a gamma ray flux higher than 
that expected by the model, indicating in particular a possible 
hardening of the spectrum during the second flare. 

1 Public quick-look results I http://xte.mit.edu/asmlc/ASM.html I 



Considering the first flare, the integral flux measured by 
ARGO-YBJ above 1 TeV is about 1.5 times higher than 
the model based on the VERITAS measurement, but still 
marginally consistent with it. The apparent disagreement be- 
tween ARGO-YBJ and VERITAS can be likely attributed to 
the non-coincidence of the data taking periods of the two de- 
tectors and to the well known small variability time scale of 
the source. The VERITAS data refer to June 6, while the 
ARGO-YBJ data are integrated over 3 days, from June 4 to 
6. Furthermore, given the difference in longitude of the two 
detectors (^160°) and the fact that they observe the source 
during few hours around the culmination time, they can never 
observe simultaneously the same object. 

The disagreement of our data with the model is more sig- 
nificant for the second flare. In order to evaluate the spectral 
behaviour in this period, we assume again a power law spec- 
trum multiplied by the EBL absorption factor e" r(£) . From our 
fitting procedure we obtain: dN/dE = (3.2±1.0) x 10 _11 (E/2.5 

TeV)- 21 ^5 e" r(£) photons cm -2 s _1 TeV -1 . 

This spectrum is shown in Fig. 3 as a solid line. The shaded 
band in the figure represents the 1 a statistical error. The sys- 
tematic errors are mainly related to the background evalua- 
tion, as discussed previously, and to the uncertainty in the ab- 
solute energy scale. According to our estimate, they globally 
affect the quoted fluxes for < 30%. 

Due to the low statistics, our data cannot constrain the 
shape of the spectrum above ^8 TeV. Nevertheless the ob- 
tained flux appears, for energies >2 TeV, significantly larger 
than predicted by Donnarumma et al. (2009), while the spec- 
trum slope is consistent with that measured by the Whipple 
Cherenkov telescope during the 2000/2001 observing season 
for a flare of comparable intensity (~ 7 times the Crab Neb- 
ula flux), also shown in Fig. 3 (dataset III, Rrennrich et al. 
(2002)). 

The integral flux measured above 1 TeV during June 11-13 
is ^6 times larger than the Crab one, making this flare one of 
the most powerful ever observed from Mrk421. 

6. DISCUSSION AND CONCLUSIONS 

Mrk421 has been continuously monitored by ARGO-YBJ 
since December 2007, showing an average VHE flux about 
twice the Crab Nebula one from February to June 2008, and 
decreasing afterwards. 

Two strong flares in June 2008 have been observed in a mul- 
tiwavelength campaign from optical to TeV energies. ARGO- 
YBJ measured the spectra of Mrk421 above 0.3 TeV dur- 
ing the second flare, completing the multifrequency obser- 
vations. A clear correlation between the gamma ray inten- 
sity measured by ARGO-YBJ and the X-ray flux measured 
by RXTE/ASM is found. 

The ARGO-YBJ data, although averaged over 3 days, ap- 
pear to support in both episodes a gamma ray flux higher 
than that predicted in the analysis of Donnarumma et al. 
(2009). However, considering the short time scale variabil- 
ity of Mrk421, it has to be noticed that our observation time 
is not fully coincident with the period referred to by the theo- 
retical curves (June 6 for the first flare and June 12-13 for the 
second one). 

The intensity of the second flare allows us to assess its spec- 
tral shape. The deabsorbed spectrum can be fitted by a power 

law oc E~ 21 -i>5 extending up to several TeV. This spectrum ap- 
pears definitively harder than that predicted on the basis of 
June 12-13 data collected up to GeV energies. 
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On the contrary, our data follow the behaviour of the energy 
spectra measured during different activity states by the Whip- 
ple Cherenkov telescope. In particular, the ARGO-YBJ data 
fully satisfy the relation between the spectral index and the 
flux obtained analyzing the measurements of Mrk421 since 
1995 (see Fig. 3 of Krennrich et al. (2002)). This result indi- 
cates that this correlation is a long term property of the source, 
as previously suggested by the Whipple collaboration. 

A global analysis of all the data collected during the 2008 
June 11-13 flare, including the present findings, could be used 
to check the compatibility of the observed phenomenology 
with current models for VHE photon emission in the jets of 



AGNs. 
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FIG. 1 . — Upper panel: rate of excess events with N pa( ; >100 observed from Mrk421 by ARGO-YBJ as a function of time from 2008 June 3 00:00 UT to June 
15 00:00 UT. Each bin contains the rate averaged over the 3 days interval centered on that bin. Lower panel: daily counting rate of RXTE/ASM. 




FIG. 2. — Rates of the excess events observed from Mrk421 by ARGO-YBJ as a function of the event minimum pad multiplicity on 2008 June 4-6 and June 
11-13 (triangles and circles, respectively). Expected rates according to the Donnarumma et al. (2009) model for June 6th and June 12-13 (dashed and solid lines, 
respectively). 
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FIG. 3. — Gamma ray flux from Mrk421 measured by ARGO-YBJ on 2008 June 11-13 (solid line). The shaded band represents 1 standard deviation error. 
The dot-dashed line sh ows the flux accordin g to the model by Donnarumma et al.(2009) for the second flare (June 12-13). The dotted line shows the spectrum 
measured by Whipple ( Krennrich et al. 2002) during a previous flare of similar intensity (see text for details). 



